Opening and exploring data
Open soil and elevation data:
#opening csv
soils <- read_csv(here("data", "shift_soil_data.csv")) %>% # read in soil csv
mutate(date = lubridate::mdy(date), # change to date class
week = lubridate::week(date)) %>% # create a week column
select(soil_id:transect, meter_location, electro_cond_mS_per_cm, season, latitude, longitude)%>% #select certain column to work with
mutate(paired_flight = lubridate::ymd(paired_flight)) %>%
drop_na() %>% #drop any NA rows
st_as_sf(coords = c("longitude", "latitude")) %>%
st_set_crs(value = 4326) %>%
st_transform(crs = 32611)
# Open spatial Data
soils_sf <- read_sf(here("data", "soil_w_elevation.shp")) %>% # read in sf file
mutate(date = lubridate::ymd(date), #change date to date class
week = lubridate::week(date)) %>% #create week
select(electro_co, date, transect, meter_loca) %>% # select certain columns
filter(date != "2022-09-12") %>%
mutate(electro_co = as.numeric(electro_co)) %>% #change to numeric class
drop_na() #drop NAs
bbox <- read_sf(here("data", "bbox.shp")) #read in bounding box file
boundary <- read_sf(here("data", "boundary.shp")) %>%
st_transform(crs = 32611)
dem <- read_stars(here("data", "dem.tif")) %>% # read in dem
st_crop(y = st_bbox(bbox)) %>% #crop to bounding box
st_warp(cellsize = 4.8, crs = 32611)
mllw <- dem + 0.042
mllw_all <- c(mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
mllw_extract <- mllw %>% # call dem
st_extract(soils$geometry) # extract elevation to point layer
Join soil and elevation data:
soils <- cbind(soils, mllw_extract$dem.tif) %>% #bind extracted values to the soil data
rename("elevation" = "mllw_extract.dem.tif") %>%
st_transform(crs = 32611)
Descriptive statistics:
#make descriptive stat table
soils %>% # to soils
group_by(transect) %>%# group by transect
summarize(min = min(electro_cond_mS_per_cm),#calculate mins
max = max(electro_cond_mS_per_cm),#calculate max
mean = mean(electro_cond_mS_per_cm),#calculate mean
dev = sd(electro_cond_mS_per_cm),#obtain standard dev
variance = var(electro_cond_mS_per_cm)) %>% # calculate variance
kableExtra::kable() %>% # create a table format
kableExtra::kable_classic(lightable_options = "striped")#theme the table
|
transect
|
min
|
max
|
mean
|
dev
|
variance
|
geometry
|
|
COPR_1_EW
|
0.42
|
1.34
|
0.6941176
|
0.2020538
|
0.0408257
|
MULTIPOINT ((235669.2 38115…
|
|
COPR_2_EW
|
6.32
|
10.53
|
8.2875000
|
1.2435749
|
1.5464786
|
MULTIPOINT ((235855 3812030…
|
|
COPR_2_NS
|
6.30
|
13.54
|
9.7581250
|
2.1347107
|
4.5569896
|
MULTIPOINT ((235855 3812030…
|
|
NCOS_1_NS
|
0.12
|
7.47
|
2.7527778
|
2.0644399
|
4.2619121
|
MULTIPOINT ((235777.4 38125…
|
|
NCOS_2_EW
|
0.74
|
8.62
|
2.9125714
|
1.8626205
|
3.4693550
|
MULTIPOINT ((235671.2 38123…
|
|
NCOS_2_NS
|
2.88
|
10.14
|
4.6538889
|
1.4925439
|
2.2276873
|
MULTIPOINT ((235685.7 38123…
|
# make descriptive stat dataframe
soil_stats <- soils %>% #call soils
group_by(transect) %>% #group by transects
summarize(mean = mean(electro_cond_mS_per_cm),# create dataframe with these stats
max = max(electro_cond_mS_per_cm),
min = min(electro_cond_mS_per_cm),
range = max - min)
Plot salinity v elevation
ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#ggplot of soils/dem data
geom_point()+#geometry of the plot
scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+# colors to use
labs(x = "Elevation (m)",#labels
y = "Electrical conductivity (mS/cm)")+
theme(legend.position = "none")+#no legend theme
guides(color = guide_legend(title = "Transect ID"))+#changing legend title
ggtitle("Soil Electrical Conductivity Across Elevations")+# title of plot
theme(plot.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),# plot theming
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 14),
axis.title = element_text(color = "#5b4f41", size = 16),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"))

Salinity vs. elevation, facet wrapped
strip_labels <- c("COPR_1_EW" = "COPR 1 EW",#create better formated labels for facet wrap strips
"COPR_2_EW" = "COPR 2 EW",
"COPR_2_NS" = "COPR 2 NS",
"NCOS_1_NS" = "NCOS 1 NS",
"NCOS_2_EW" = "NCOS 2 EW",
"NCOS_2_NS" = "NCOS 2 NS")
ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#start ggplot
geom_point()+# point geometry
scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+#colors
facet_wrap(~transect, labeller = as_labeller(strip_labels))+#subplots based on transect
theme(legend.position = "none")+# legend theme
labs(x = "Elevation (m)",#labels
y = "Electrical conductivity (mS/cm)")+
guides(color = guide_legend(title = "Transect ID"))+#legend title
ggtitle("Soil Electrical Conductivity Across Elevations")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot theme
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"))

Open Raster Files:
Opening and Compiling NDVI:
ndvi_1 <- read_stars(here("imagery", "NDVI", "2022_02_24_ndvi.tif")) %>%
setNames("ndvi")
ndvi_2 <- read_stars(here("imagery", "NDVI", "2022_02_28_ndvi.tif")) %>%
setNames("ndvi")
ndvi_3 <- read_stars(here("imagery", "NDVI", "2022_03_08_ndvi.tif")) %>%
setNames("ndvi")
ndvi_4 <- read_stars(here("imagery", "NDVI", "2022_03_16_ndvi.tif")) %>%
setNames("ndvi")
ndvi_5 <- read_stars(here("imagery", "NDVI", "2022_03_22_ndvi.tif")) %>%
setNames("ndvi")
ndvi_6 <- read_stars(here("imagery", "NDVI", "2022_04_05_ndvi.tif")) %>%
setNames("ndvi")
ndvi_7 <- read_stars(here("imagery", "NDVI", "2022_04_12_ndvi.tif")) %>%
setNames("ndvi")
ndvi_8 <- read_stars(here("imagery", "NDVI", "2022_04_20_ndvi.tif")) %>%
setNames("ndvi")
ndvi_9 <- read_stars(here("imagery", "NDVI", "2022_04_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_10 <- read_stars(here("imagery", "NDVI", "2022_05_03_ndvi.tif")) %>%
setNames("ndvi")
ndvi_11 <- read_stars(here("imagery", "NDVI", "2022_05_11_ndvi.tif")) %>%
setNames("ndvi")
#ndvi_12 <- read_stars(here("imagery", "NDVI", "2022_05_12_ndvi.tif")) %>%
# setNames("ndvi")
ndvi_13 <- read_stars(here("imagery", "NDVI", "2022_05_17_ndvi.tif")) %>%
setNames("ndvi")
ndvi_14 <- read_stars(here("imagery", "NDVI", "2022_05_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_all <- c(ndvi_1, ndvi_2, ndvi_3, ndvi_4, ndvi_5, ndvi_6, ndvi_7, ndvi_8, ndvi_9, ndvi_10, ndvi_11, ndvi_13, ndvi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
Opening and Compiling VSSI:
vssi_1 <- read_stars(here("imagery", "VSSI", "02_24_vssi.tif")) %>%
setNames("vssi")
vssi_2 <- read_stars(here("imagery", "VSSI", "02_28_vssi.tif")) %>%
setNames("vssi")
vssi_3 <- read_stars(here("imagery", "VSSI", "03_08_vssi.tif")) %>%
setNames("vssi")
vssi_4 <- read_stars(here("imagery", "VSSI", "03_16_vssi.tif")) %>%
setNames("vssi")
vssi_5 <- read_stars(here("imagery", "VSSI", "03_22_vssi.tif")) %>%
setNames("vssi")
vssi_6 <- read_stars(here("imagery", "VSSI", "04_05_vssi.tif")) %>%
setNames("vssi")
vssi_7 <- read_stars(here("imagery", "VSSI", "04_12_vssi.tif")) %>%
setNames("vssi")
vssi_8 <- read_stars(here("imagery", "VSSI", "04_20_vssi.tif")) %>%
setNames("vssi")
vssi_9 <- read_stars(here("imagery", "VSSI", "04_29_vssi.tif")) %>%
setNames("vssi")
vssi_10 <- read_stars(here("imagery", "VSSI", "05_03_vssi.tif")) %>%
setNames("vssi")
vssi_11 <- read_stars(here("imagery", "VSSI", "05_11_vssi.tif")) %>%
setNames("vssi")
#vssi_12 <- read_stars(here("imagery", "VSSI", "05_12_vssi.tif")) %>%
# setNames("vssi")
vssi_13 <- read_stars(here("imagery", "VSSI", "05_17_vssi.tif")) %>%
setNames("vssi")
vssi_14 <- read_stars(here("imagery", "VSSI", "05_29_vssi")) %>%
setNames("vssi")
vssi_all <- c(vssi_1, vssi_2, vssi_3, vssi_4, vssi_5, vssi_6, vssi_7, vssi_8, vssi_9, vssi_10, vssi_11, vssi_13, vssi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
Opening and Compiling mARI:
mari_1 <- read_stars(here("imagery", "mARI", "02_24_mari_avg.tif")) %>%
setNames("mari")
mari_2 <- read_stars(here("imagery", "mARI", "02_28_mari_avg.tif")) %>%
setNames("mari")
mari_3 <- read_stars(here("imagery", "mARI", "03_08_mari_avg.tif")) %>%
setNames("mari")
mari_4 <- read_stars(here("imagery", "mARI", "03_16_mari_avg.tif")) %>%
setNames("mari")
mari_5 <- read_stars(here("imagery", "mARI", "03_22_mari_avg.tif")) %>%
setNames("mari")
mari_6 <- read_stars(here("imagery", "mARI", "04_05_mari_avg.tif")) %>%
setNames("mari")
mari_7 <- read_stars(here("imagery", "mARI", "04_12_mari_avg.tif")) %>%
setNames("mari")
mari_8 <- read_stars(here("imagery", "mARI", "04_20_mari_avg.tif")) %>%
setNames("mari")
mari_9 <- read_stars(here("imagery", "mARI", "04_29_mari_avg.tif")) %>%
setNames("mari")
mari_10 <- read_stars(here("imagery", "mARI", "05_03_mari_avg.tif")) %>%
setNames("mari")
mari_11 <- read_stars(here("imagery", "mARI", "05_11_mari_avg.tif")) %>%
setNames("mari")
#mari_12 <- read_stars(here("imagery", "mARI", "05_12_mari_avg.tif")) %>%
# setNames("mari")
mari_13 <- read_stars(here("imagery", "mARI", "05_17_mari_avg.tif")) %>%
setNames("mari")
mari_14 <- read_stars(here("imagery", "mARI", "05_29_mari_avg.tif")) %>%
setNames("mari")
mari_all <- c(mari_1, mari_2, mari_3, mari_4, mari_5, mari_6, mari_7, mari_8, mari_9, mari_10, mari_11, mari_13, mari_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
Opening and Compiling CRSI:
crsi_1 <- read_stars(here("imagery", "CRSI", "02_24_crsi.tif")) %>%
setNames("crsi")
crsi_2 <- read_stars(here("imagery", "CRSI", "02_28_crsi.tif")) %>%
setNames("crsi")
crsi_3 <- read_stars(here("imagery", "CRSI", "03_08_crsi.tif")) %>%
setNames("crsi")
crsi_4 <- read_stars(here("imagery", "CRSI", "03_16_crsi.tif")) %>%
setNames("crsi")
crsi_5 <- read_stars(here("imagery", "CRSI", "03_22_crsi.tif")) %>%
setNames("crsi")
crsi_6 <- read_stars(here("imagery", "CRSI", "04_05_crsi.tif")) %>%
setNames("crsi")
crsi_7 <- read_stars(here("imagery", "CRSI", "04_12_crsi.tif")) %>%
setNames("crsi")
crsi_8 <- read_stars(here("imagery", "CRSI", "04_20_crsi.tif")) %>%
setNames("crsi")
crsi_9 <- read_stars(here("imagery", "CRSI", "04_29_crsi.tif")) %>%
setNames("crsi")
crsi_10 <- read_stars(here("imagery", "CRSI", "05_03_crsi.tif")) %>%
setNames("crsi")
crsi_11 <- read_stars(here("imagery", "CRSI", "05_11_crsi.tif")) %>%
setNames("crsi")
#crsi_12 <- read_stars(here("imagery", "CRSI", "05_12_crsi.tif")) %>%
# setNames("crsi")
crsi_13 <- read_stars(here("imagery", "CRSI", "05_17_crsi.tif")) %>%
setNames("crsi")
crsi_14 <- read_stars(here("imagery", "CRSI", "05_29_crsi.tif")) %>%
setNames("crsi")
crsi_all <- c(crsi_1, crsi_2, crsi_3, crsi_4, crsi_5, crsi_6, crsi_7, crsi_8, crsi_9, crsi_10, crsi_11, crsi_13, crsi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
Compiling NIR Bands
nir_1 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_24.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_2 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_28.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_3 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_08.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_4 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_16.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_5 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_22.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_6 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_05.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_7 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_12.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_8 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_20.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_9 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_10 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_03.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_11 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_11.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
#nir_12 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_12.tif")) %>%
# st_set_dimensions(3, values = seq(1:425)) %>%
# filter(band %in% c(76:126)) %>%
# split(3) %>%
# setNames(c(76:126))
nir_13 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_17.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_14 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
filter(band %in% c(76:126)) %>%
split(3) %>%
setNames(c(76:126))
nir_all <- c(nir_1, nir_2, nir_3, nir_4, nir_5, nir_6, nir_7, nir_8, nir_9, nir_10, nir_11, nir_13, nir_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date") %>%
st_warp(dest = crsi_all)
Adjusting elevation data
mllw_all <- mllw_all %>%
st_warp(dest = crsi_all) %>%
setNames("elevation")
PCA
PCA Object
# Read in the data
nir_pcd <- soils_extract %>% # call object
as.data.frame() %>%
select(electro_cond_mS_per_cm,
`76`:`126`) %>% # select the variables we want to explore
drop_na() # drops columns with NAs as PCA only works with numeric values
Generate Biplot
#Creating the biplot via autoplot()
autoplot(nir_pca, # pca data
data = nir_pcd, # original data
loadings = TRUE, # show loadings
loadings.label = TRUE, # label loadings
loadings.colour = "black", # loading color
loadings.label.colour = "black", # loading label color
loadings.label.vjust = -.5) + # vertical justification
ggtitle("Principle Component 1 and 2 Biplot (SILVA)") + # title
labs(x = "Principle Component 1",
y= "Principle Component 2")+ # axis labels
theme(plot.title = element_text(color = "#5b4f41", size = 16),
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash",
color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("white"),
axis.line = element_line(color = "#5b4f41")) # change themes

# Creating screeplot
sd_vec <- nir_pca$sdev # standard deviation
var_vec <- sd_vec^2 # variance
pc_names <- colnames(nir_pca$rotation) # names of PCs
pct_expl_df <- data.frame(v = var_vec, # new data frame for screeplot; variance
pct_v = var_vec / sum(var_vec), # percent of variance
pc = fct_inorder(pc_names)) %>% # orders rows
mutate(pct_label = paste0(round(pct_v * 100, 1), "%")) # adds character %
ggplot(pct_expl_df, aes(x = pc, y = v))+ # graphs via PC and variance
geom_col()+
geom_text(aes(label = pct_label), vjust = 0, nudge_y = 0.008)+
labs(x = "Principle Components",
y = "Variances") +
ggtitle("PCA Screeplot")+
theme(plot.title = element_text(hjust = .5))

Fitting model with all variables currently available
Fitting the model to all variables
fit_all <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi,
data = soils_extract,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
Evaluating Model
# using caret
fit_all_imp <- as.data.frame(fit_all$finalModel$importance)
fit_all_imp_scaled <- predict(preProcess(fit_all_imp, method = c("range")), fit_all_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_all_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

cor(fit_all$finalModel$y, fit_all$finalModel$predicted)
## [1] 0.8421477
Ploting Model Fit
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_all$finalModel$y, x = fit_all$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation and Indices (caret method)")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted)
all_corr
## [1] 0.8421477
Fitting the model to spectral Indices
fit_si <- train(electro_cond_mS_per_cm ~ ndvi + mari + vssi + crsi + `76` + `77` + `78` + `79` + `80` + `81` + `82` + `83` + `84` + `85` + `86` + `87` + `88` + `89` + `90` + `91` + `92` + `93` + `94` + `95` + `96` + `97` + `98` + `99` + `100` + `101` + `102` + `103` + `104` + `105` + `106` + `107` + `108` + `109` + `110` + `111` + `112` + `113` + `114` + `115` + `116` + `117` + `118` + `119` + `120` + `121` + `122` + `123` + `124` + `125` + `126`,
data = soils_extract,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
Evaluating Model
fit_si_imp <- as.data.frame(fit_si$finalModel$importance)
fit_si_imp_scaled <- predict(preProcess(fit_si_imp, method = c("range")), fit_si_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_si_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_si_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_si_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

cor(y = fit_si$finalModel$y, x= fit_si$finalModel$predicted)
## [1] 0.5060143
Observed v Predicted Plot
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Spectral Indices")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "si_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
Random Forest regression of Soil salinity and elevation
fit_ele <- train(electro_cond_mS_per_cm ~ elevation, data = soils, method = "rf", trControl = trainControl(method = "cv", number = 10), importance = T)#random forest fitting of the data with stated equation
Correlation
elevation_cor <- cor(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted)#calculate correlation coefficient
elevation_cor#call/print correlation value
## [1] 0.9208296
Observed v Predicted
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "ele_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
Making a Map
all_pred_feb <- all_pred %>%
filter(date == "2022-02-24") %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE)
ggplot()+
geom_stars(data = all_pred_feb, aes(x, y, fill = prediction)) +
scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
labs(y = "Northing",#labels
x = "Easting",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("Soil Salinity")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))

TESTING OF PLSR INSTEAD!
Fitting the model to all variables
fit_all_pls <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi,
data = soils_extract,
method = "pls",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
Evaluating Model
# using caret
predicted_values <- as.data.frame(fit_all_pls$finalModel$fitted.values) %>%
select(`.outcome.3 comps`)
training_data <- as.data.frame(fit_all_pls$trainingData$.outcome)
train_pred <- c(predicted_values, training_data) %>%
as.data.frame()
cor(training_data, predicted_values)
## .outcome.3 comps
## fit_all_pls$trainingData$.outcome 0.6955281
Ploting Model Fit
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color = "#3A5D3D")+#smooth line data
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation and Indices (PLSR)")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_pred_pls <- predict(raster_all, fit_all_pls, drop_dimensions = F)
all_pred_feb_pls <- all_pred_pls %>%
filter(date == "2022-02-24") %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE)
ggplot()+
geom_stars(data = all_pred_feb_pls, aes(x, y, fill = prediction)) +
scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
labs(y = "Northing",#labels
x = "Easting",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("Soil Salinity")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
